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NOMENCLATURE 


F(5) 

H(x) 

hi 

K 

M 

M 

no 

n 

p 

Pi 

Q 


thickness distribution 
distribution of camber 

half of the height of 2h^ = rj + n i 

1 - 

transonic similarity parameter, K = 

M 

numbe-r of intervals into which tho airfoil is divided by 
rectangular elements 

free- stream Mach number 

number of rectangular elements or mesh points 

- (x.y) 

= 

- u,o 

magnitude of the difference between yj^ and the y-coordinate 
of the bottom edge of Aj 

magnitude of the difference between yj^ and^ tiie y-coordinate 
of the top edge of Aj^ 


II, V 

u, V 

‘^i 

x,y 

x.y 

Xf.yi 

Y^(x), Y 

Y^(x) 


velocity components in the tangential and transverse directions, 
— 3d) — 3d) 

” = , V = -ri 

3x 3y 


respectively, u 


2$ 


velocity components in the transformed plane, ii - , 

tangential velocity at the 1th mesh point 

rectangular Cartesian coordinates in the physical plane 

rectangular Cartesian coordinates in the transformed plane, 

X = X, y =* By 

coordinates of the 1th mesh point 
(x) upper and lower airfoil profiles, respectively 


3 ^ 
3y 


modified airfoil profile, Y^(x) 


ili 


Y + 1 1 - 
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e 
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U 

T 

^(x,y) 

<l>(x,y) 

^L^x.y) 

♦Lj 

^jj(x.y) 


nngU' of attack, deg 

Fraud t] -Clanert parameter, P = 

ratio of specific heats 

vortex strength, y(x) = ti(x,+0) - u(x,-0) 
= Y(xi) 

= u^(x,y) - u^(x,-y) 
dY^ dY_ 

" "dT ■ 


the ith rectangular element 
thickness ratio 
half the width of 
half the width of a 

e 

aspect ratio of 

. 2(y -f 1) 


rectangular Cartesian coordinates 
the vanishing rectangular cavity 
measure of the distribution of camber 
perturbation velocity potential 
modified perturbation velocity potential. 


♦(x,y) = 


■(y + l)K.n_ _ 
_ — 


y) 


contribution to potential due to lifting effects 


= <l't (x^.y^) 


contribution to potential due to nonlifting effects 
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SOLUTION OF TRiXNSONIC FLOWS BY AN 


INTEGRO-DIFFEKFNTIAL EQUATION METHOD 
Wandera Ogana* 

Ames Research Center 


SUMMARY 


For free-stream Mach numbers less than unity, steady transonic flow past 
a two-dimensional airfoil is described by a singular integro-dif f erentia 1 
equation which Involves the tangential derivative of the perturbation veloc- 
ity potential. The integro-dif ferential equation can be transformed into a 
nonlinear system of algebraic equations involving numerical evaluation of a 
first-order derivative. Subcrltical flows are solved by a central difference 
approxinuition of the derivative everywhere. For supercritical flows with 
shocks, central differences are taken in locally subsonic flow regions and 
backward differences in locally supersonic flow regions. A normal shock con- 
dition is applied. A condition is enforced at the sonic point which ensures 
smooth transition from subsonic to supersonic velocities. The region of Inte- 
gration is partitioned into rectangular elements whose sizes can be adjusted 
in order to minimize the number of mesh points. The method is applied to a 
nonlifting parabolic-arc airfoil in both subcrltical and supercritical flows, 
and to a subcrltical flow over a lifting NACA 0012 airfoil. Although a small 
number of mesh points is used, the results compare favorably with finite dif- 
ference results. For lifting flown, a modified boundary condition should be 
applied to Improve the accuracy of the results at the leading edge of the 
airfoil . 


INTRODUCTION 


The transonic small perturbation equation for steady, two-dimensional 
flows is transformable to a singular integro-dif ferential equation which can 
be differentiated to yield a singular integral equation in the tangential 
velocity (refs. 1 and 2). Solutions of the Integral equation have been 
obtained for subcrltical flows (refs. 3-6) and for embedded supersonic regions 
with shock waves (refs. 3, 4 and 7). The transonic integro-dif ferential equa- 
tion has not been solved numerically. A related linear integral equation for 
three-dimensional subsonic flows has been studied by Morino and Kuo (ref. 8). 
The present method solves the integro-dif ferential equation using a technique 
which Ogana (refs. 6 and 9) developed for the transonic integral equation. 
Solutions for supercritical flows with shocks can be determined by taking 
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lontriil H i f f ereiu'L's in local subsonic flow regions and backward differemes 
in local supersonic flow regions as propi)sed by Murman and Cole (ref. 10). 
Special considerar ion is given to sonic and shock points. An advantage of 
the int egro -di f f erent ia i equation metliod over the integral equation method is 
that the shock conditions can be introduced in a simple direct mam. 

INTEGRO-DIFFERENTIAL EQUATION 

We consider steady, two-dimensional transonic flow past a thin airfoil 
of thickness ratio 6, and measure of the amount of camber t, at a small 
angle of attack o. For subsequent convenience of notation, the rectangular 
Cartesia.i coordinates in the physical plane are taken to be (x,y). Let the 
upper airfoil profile be described by 

Y^(x) * 6F(x) + tH(x) - ax (la) 

and tlie lower profile by 

Y_(x) = -6F(x) + tH(x) - ax (lb) 

where F(x) and H(x) are the thickness and camber distributions, respectively. 

Let ^(x,y) be the perturbation velocity potential. The perturbation 
velocities in the tangential and transverse directions are given by u = i- 
and V » respectively. The transonic small perturbation equation for 
^ is 

[1 - M ^ - (y + DM 2 J = 0 (2) 

oo ■' 00 XX yy 

where is the free-streain Mach number and y is tne ratio of specific 

heats. 

Equation (2) is solved subject to the Kutta condition, the body boundary 
condi tlon , 

dY* 

v(x,±0) = (3) 

and the condition that u and v vanish at infinity. 

The following shock jump conditions hold (ref. 10): 

M + D 

o._ ® — ^ ... 


<(1 - M„2)G - 

— G^) (dy)j, - 

( V ) (dx)j. - 0 

(4a) 


< 5> (dy)j. + 

( u ) (dx) j. - 0 

(4b) 


where ( ) denotes the jump in a quantity across the shock, and the subscript 
E refers to an element on the shock surface. 
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Th 

e transformations, » 

1 - M.„2, 

X 
II 

XI 

y “ 

f’y. 


>Hx,v) 

= I(> + 1)M„,^/B’’]5(x,y), 

u(x,y) * 

•I't /^X 

and 

v(x,y) = a4>/3y. 


reduce 

(eqs. (2) .md (3)) to 







(1 - 

'<>x)''<'xx 


0 


(5) 










v(x 

,t0) = 

dx 


(6) 


wliere Y+(x) * ( (y + 1 ) /6K ^ ] Y+ (x) defines the modified airfoil profile and 

K = (1 - ^ is the transonic similarity parameter. 

For a normal shock, equations (4a) and (4b) become 


(u--|^u^)>"0 (7) 


Let (K,K) he the running coordinates in the integration; then, for 
< 1, application of Green's theorem converts (eqs. Cj)-(7)) to the follow- 
ing singular Integro-dif ferential equation (refs. 1, 2 and 11): 


())(x,y) = 4>j^(x,y) + (}>j^(x,y) 



(C.OdS 


( 8 ) 


where and -/t are defined below. 


The contribution to the potential due to nonlifting effects is 


<tfj(x,y) 




- x)2 + y2] 


1 /. 


d^ 


where 


A4.^(C) 


dY^ dY 

’T 

ir * "di 


(9) 


( 10 ) 


The contribution to the potential due to lifting effects is 
i-an be written in the following form after an Integration by parts (ref. 11): 

4>L(x,y) " jfsgn(y) - arc tan j Y (Od^ (11) 


where 

yiO - u(C,+0) - u(C.-O) (12) 


defines the local strength of the vortices distributed along the chord. For 
a nonlifting airfoil, “ 0, hence <^j^(x,y) ■ 0. 
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Thf kernel of the Integro-di f f erentlal equation Is 


.X' (C - X, ; - y) 


- X)? +V^’ - y)-'') 


( 13 ) 


The singular point of the kernel is enclosed in a vanishing rectangular 
cavity. Op, of aspect ratio X (fig. 1). The significance of this cavity was 
reported in a previous article by Ogana and Spreiter (ref. 2). 


NONLIFTING EFFECTS 


From equations (1) and (10) we deduce that A(>^((;) ■ ( 2(\ + 1)/K^^^ )F' (^) . 
If this relation is substituted in equation (9) and an integration by parts 
performed, then for an airfoil of unit chord we obtain 




Y 1 X - ^ 

(C - x)2 + y2 


F(5)dC 


(14) 


LIFTING EFFECTS 


For lifting airfoils, equation (8) can be solved only after y( 0 in 
equation (11) has been determined. To do this we differentiate both sides of 
equation (8) with respect to y to find 


v(x 


.^) • ^ r 




(C - X)2 + y2 C 
Xi'(c - X, c - y) 


HriOdi + it 


t/' 


C - X 


[// 


2 Vo (^ - X)2 + y2 

u2 (t,Odsj 


Y(OdC 


(15) 


We now take the limit as y -*• +0 and y -0 on both sides of equation 
(15), add the two results and perform an integral inversion. The procedure, 
described in appendix A, yields the result 


where 


y(x) 




u - 2 (y + 


constant 


( 16 ) 


4 


o •> o 

Au‘(x,y) * u (x,y) - n (x, -y) 


'•(O = X^x- C,y)Au^(x,y)dx dy 

S 

From equation (A6) we can deduce 


X(x - ^,y) 


y(C - x) 


4it[(^ - x)2 + y2)2 


(17) 

(18) 


(19) 


MATRIX EQUATION 


Let P = (x,y), Q = (C,C) and .lli'v'Q.P) = - x, r, - y). Equation 

(8) can be written 

*(P) = «l»j^(P) + yyJi^(Q.P)u^(Q)ds 


( 20 ) 


where 


u(Q) - -^(♦(Q)] 


( 21 ) 


The region of integration, S, is now divided into n rectangular ele- 
ments, labelled Aj^, A 2 , . . ., Aj^. In each Aj there is a mesh point 
Pi(xi,yi), located as shown In figure 1. The mesh points are concentrated 
near the airfoil and become fewer as the elements grow larger (fig. 2). It 
is now assumed that the flow quantities are constant within each element. 
For a mesh point, Pj, equation (20) becomes 


'(Pi> = <l>j,(Pi) + 4>L(Pi) +2 y!xr (Q,Pj)u2(Q)dS 


( 22 ) 


.1*^1 A 


J 


Unless otherwise stated it will be assumed that the subscripts i, j * 1, 
2, . . ., n. 


Let 


= 4 »(Pj). '(•Ni = ^ 


(23) 


and define *:he matrix (b^j ] by 


'ij 


yjir(Q,p^ 


)dS 


(24) 
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Tlien equation (22) beeomes 




i 


<f>V + <}'I . 

•’i *1 



( 25 ) 


When the right-hand side of equation (21) is ev.jluated numerically, U| 
is expressed as a function of the discrete values of the perturbation poten- 
tial. Equation (25) can then be regarded as a system of nonlinear algebraic- 
equations in ^2* • • • ♦ <^n • Dt^fine the component gj(^) of the vector 

g(T) by 


8j(v) 



(26) 


In matrix form, equation (25) becomes 

I * g(t) (27) 

where 

♦ * <l>2» • • •* (28a) 

g(l) = (gi(^). g^(|). • . •. gn(<>)l^ (“8b) 


Equation (14) shows that for a given airfoil protiie, can be evalu- 

ated immediately by some suitable numerical integration scheme. Before solu- 
tions of equation (27) can be obtained, certain matrices have to be evaluated. 
Suppose the airfoil is divided into M Intervals by the rectangular elements. 
Using equation (11) we can write as 

M 

k-1 

where 


Yk - Y(x^) 


’ik 


?.v J . 2 


: 5 -sgn(yj^) - arc tan 



(30) 


The matrix elements, q^j^t ar** evaluated in sec. I of appendix B. 

To obtain we let 12(x,y) - Au^(x,y) and *- h(xj^,y|^). Equation 

(16) yields 
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(31 ) 


where 






( i '2 ) 


Expressions for the matrix elements, derived in see. 2.1 of 

appendix B. 

In order to evaluate in equation (31), we require Gi, (J?, . . ., 

Using equations (18) and (19) we can obtain 


where 


fj 


// 


n 

j = l 


y(x^ - x) 

4tt[ (Xf - x)^ 


(33) 


dx dv 


(3i) 


The matrix elements, d^^^, are determined in see. 2.2 of appendix B. 

Tb® matrix defined by equation (24) is associated with the nonlinear eon- 
tributton to the potential. It can be written as 


b. . = 
ij 


// 


Xi - ? 


'in [ (C - Xj^)^ + (4 - yi)“] 


^ dC dC 


(3 3) 


dix B. 


Expressions for the matrix elements, b^j , arc derived in sec. 3 of appen- 


SOLUTION OF THE MATRIX EQUATION 


Define <<>n ’ * ' *’ Equation (27) is solved by the 

Jacobi iteration as follows: 


<P 

>(m+l ) 


N 


m = 0, 1, 2, . . 

where the superscript refers to the iteration number. 


(36a) 

(36b) 
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K'lu*n siiIvLnp oqu.it ion (27) we .ire eoni' ron ted witli 


II . 
1 



■ iiul 


07 ) 


Numerical evaluation of this derivative depe.nds on tlie nature I’f tlie 
It'cal liow as determined by the sipn of tlie coefficient of in equation 

(S). When (1 - 4^) ' 0, equatii'n (S) Is elliptic and describes a local sub- 
sonic flow; wlien (1 - i>y) = 0, it is parabolic and ilescribes a local si'iiic 
flow; and when (1 - 4>^) • 0, it is livperbolic and describes a loc.il supersonic 
flow. The difference scheme proposed by Murman and i\'le (ref. 10) is designed 
to solve tin’s tvpe of mixed flow equation. Phe main feature of the scheme is 
th.it, in the tangential di-ection, leutral differences iro taken in locally 
subsonic flow regii'ns and backward differences in locally supersonic flow 
ri'gions. We extend these ideas to the present method. 


To simplify the analysis, we assume that the r-,ctangular elements are 
■qua] within .» given strip parallel to the x-axis (fig. 2). They can differ 
in si7.e from one strip to another. The mesh points are numbered so that 
xs^j x^. Lit Ox width of each element, then Xj+I ” * ^x» 

.all X j in the strip. 


Wo define the following two operators at a mesh point Pj : 




- 4> 


izi 


2Ax 


(38) 


4>i 

V = 1 - — — 

H 2Ax 




(39) 


During an iteration, Vj? and are computed at each mesh point, If 

\’e ^ 0, tho flow at Pj is subsonic and equation (37) is evaluated by the 
central difference formula; i.e.. 


. Vi - *J-i 
j 


(40) 


If Vg < 0 and < 0, the flow at Pj is supersonic and equation (37) 
is evaluated by the backward difference formula; i.e,. 


u . 
1 


— f 


j ~ *j-2 
2Ax 


(41) 


As the flow accelerates through sonic velocity from subsonic to super- 
sonic velocities, a point is reached where Vg < 0 and > 0. Using the fact 
that for a locally sonic flow, 1 - “ 0, we set 

u. - 1 (42) 
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This is analogous Co the parabolic point operat»>r (r*. 12 and 13). 


As the flow decelerates across the shoik iron supersunii in siihsonic 
velocities, a point is reached where Vg > 0 and V|^ < 0. Mnrm.in (ri>f. 1 'O 
introduced a shock-point operator t(» be applied at such a point. We use a 
different approach since it is impractical to extend the shock-point operator 
concept to the present system. We assume that the mesh point, Pj, where 
Vg > 0 and < 0, is located just downstream of the shock. The adjac-enl 
mesh point, located upstream of the shock, is ^j-l* Prom equation (7) it 
can be deduced that 

Uj_^ + Uj = 2 (4 i) 

If 'Jj-1 is given, then uj can he determined from this equation. The 
use of Uj-1 to determine uj is based on the fact that information propa- 
gates downstream in a locally supersonic flow region. 

In equations (38) through (41), we could use unevenly spaced grids where 
appropriate and apply the difference formulas given by Bailey (ref. 15), for 
instance . 

For mesh points adjacent to the far field boundary, equation (37) is 
evaluated using a method which Ames (ref. 16) described for derivatives near 
a boundary. If the f<.. field boundary is infinitely far from the airfoil, 
the potential on it could be held fixed through all iterations. The computa- 
tional region does not, in practice, extend very far. Consequently, it is 
best to update the potential on the far field boundary after a few iterations. 

Before solving equation (27), it is important to give careful considera- 
tion to the number of mesh points required. Figure 2 shows the elements 
growing larger the farther they are from the airfoil. This type o-f partition- 
ing enables computation to be performed in a large region of integration 
while using a small number of mesh points. It also ensures that the '(Signif- 
icant contribution from mesh points near the airfoil is not dispersed in 
unusu.illy large elements. For a given distribution of the elements, tiie 
computational region can be extended by increasing the aspect ratios of the 
elements. 


RESULTS AND DISCUSSION 


Figures 3 and 4 show the coefficient of /“’essure plots for a nonlifting 
parabolic-arc airfoil of thickness ratio 6 - 0.06 at free-stream Mach num- 
bers 0.825 and 0.87, respectively. Figure 5 shows the coefficient of pressure 
plots for a NACA 0012 at a * 2® and = 0.63. The results for the present 
method are compared with those of a finite difference scheme by Ballhaus, 
Jameson and Albert (ref. 17) in figure 5. The finite difference method uses 
a 67x41 mesh point network; convergence is achieved in 40 to 100 iterations, 
depending on the free-stream Mach number. 
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Tm compnrison, the present method uses uhout 200 mesli piunts located in 
rectangular elements whose aspect ratios range from 2 to S. Iteration pro- 
ici-ds until the largest difference between successive iterates is smaller 
than some specified tolerance; i.e.. 


max 

i 


I ^ ("1+1 ) 



(44) 


where t is the tolerance. 


The rate of convergence for the present method, applied to the parabolic- 
arc airfoil, is shown in the following tabulation for two tolerance values. 


M 

00 

Tolerance 

I terations 

0.825 

10"3 

4 


10-'* 

7 

0.870 

10'3 

12 


lO"*" 

20 


The results presented are obtained when the potential on the far field 
boundary is updated after each iteration. If the potential on the boundary 
is held fixed, there is a noticeable effect on values near the boundary, but 
that effect diminishes in the interior. 

Near the leading edge of a lifting airfoil, the accuracy of the method 
could be improved by applying a modified boundary condition described by 
Nixon (refs. 7 and 18). 


CONCLUDING REMARKS 


The present study indicates that the integro-dif f erent ial equation method 
can provide rapid and reliable solutions for transonic flow problems. Con- 
trary to what happens in the integral equation schemes, the shock conditions 
are enforced in a simple and straightforward manner. Although a relativelv 
small number of mesh points is used, the results agree well with finite dif- 
ference computations. 

Extension to three-dimensional lifting flows can be made provided it is 
possible to determine the equivalent of the integral inversion that yields an 
expression for the vortex strength. There should be no fundamental problems 
in applying the present scheme to nonlifting three-dimensional flows. 
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APPENDIX A 


VORTEX STRENGTH 


The vortex strength, y(x), 
integral equation 


is determined by solving the 


following 


1 inear 


v(x,y) 







C - X 


x)‘ + v" 


v(OdC 



; - y)u"^ (C,;)dS 




(15) 


where A4>j.(^) is given by equation (10) and .>l'(C - x. ^ - y) by equation (13). 

We now take the limit as y ->■ ±0 on both sides of equation (15) and 
add the two results. 

Karamcheti (ref. 19) shows that 


lim f- 

- 10 


U - x)^ + y2 


A4)j.(4)dC * ±^A(^^(x) 


(Al) 


Mangier (ref. 20) shows that 
1 


lim ^ I ■ 
r ■> ±0 


- r Y(OdC , 

u - xr + y 




y(0 

C - X 


dC 


(A2) 


where the right-hand side is a Cauchy integral. 

Nixon and Hancock (ref. 5) have evaluated the last term in equation (15) 
as y ±0. From equations (1) and (6) we deduce 


v(x,40) + v(x,-0) * u[tH'(x) - a] (A3) 

where the constant u * 2 (y + 1)/6K^^^. 

Using equation (A3) together with the results obtained when y -► lO 
on the right-hand side of equation (15) we conclude 
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f 

( 


J"L'" - 'J - - J -J-SV J T 


//■ 


+ f I - x,f.) Au^ (C, .) ds 
S 


will' re 


Au‘ {E,, O = K, C) - {E,,-C.) 


\{E - X, ;) = 


'.(x - O 


4n 


[l5 - x)^ + cO’ 


When the terms are suitably arranged, equation (A^i) can be inverted 
to give an explicit expression for >(x). We define 


AW(C) = >(fj - 


Au‘ (f , 0) 


c;(x) = //« - X, 4) Au^ (E, 4) dc dc 
S 

W^,(x) = -|ij|^rH'(x) - aj - G(x) 


Substitution of equations (.A7) - (A9) into equ t ion (A4) yields 

1 




Ashley and Landahl (ref. 21) show that equation (AlO) can be 
inverted to give 


1 


iW(x) f 


2 rr^ X / "o'"' ^-p: 


w iO 


We consider an airfoil without camber and use the following result 
from Ashley and Landahl (ref. 21): 


(A4) 

(A5) 

(A6) 

(A7) 

(A8) 

(A9) 

(AlO) 

(All) 
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IT I X - C T 1 - ^ 


-1 


Substituting for AW(x) and W_(C), we obtain the expression 




/X ^ - X , Au- (x.O) , 

^(x) = paV — - — + — + 




(A12) 


(A13) 
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APPENDIX B 


MATRIX ELEMENTS 


In tills ;ippi*ndix we derive expressions for sone matrix elements which 
are required in the numerical solution of the inte>»rn-dif ferent ial equn»'ion. 

1 . Ma trix Associated I’ith the Ton tribution of Lift to tlie Po tent ial 

The matrix associated with the contribution of lift to the potential 
at a mesh point P^Cx^.y^) is given by 


'’ik 


"k L 


sgn(y^) - 


arc tan 



(30) 


where i •-« 1,2, . . . ,n, k = 1,2, . . . ,M, and y^ 0 . 
Equation (30) can be evaluated at once to give 


^ik 2 tt 


n6j^sgn(yj^) - (Xj^ “ ’'i tan^ 


X, - X 




+ (Xj^ ~ ~ <5j^) arc tan 


Vi 

+ ^ «n 

2 


Vl -k ‘“k - "l 'k 

* <’‘k • 


(Bl) 


for y. = iO. 
1 


We have to consider separately three values of xj, namely, xj < Xj^, 
Xj^ = X|^ and x^^ > xj^. Integration of equation (30) yields the result 


‘’ik = 


0 , 


^ ’'k 


- 6 ^ sgn(y^), x^ - Xj^ 

^ 6^ sgn(yj^). 


^ \ 


(B2) 


lA 





2 . Matrices Associated with the Vortex Strength 


There are two matrices associated with the vortex strength; one arises 
from a double integral before an inversion is performed and the other arises 
from a line integral after the inversion. 


2.1 Matrix from the Line Integral 
The matrix element is given by 


C-x,Vl - C • 


'k£ 


. ,M 


(32) 


C - Xj^, b - x,^, a - 1 - x^, t 2 - - Xk + tj « x^ - x,,^ - 6^. 


Equation (32) is transformed to 


'kH 


2 /* dt 1 /b t 

IT y t K a - t 


(B3) 




The substitution r^ * (b + t)/(a - t) converts the integral to one 
which can be evaluated easily to give 



where 


T(t) 


arctan 




b + t 
a - t 




/a(b t) - /b(a - t ) 
/a(b + t) + /b(a - t) 


(B4) 


(B5) 


2 . 2 Matrix from the Double Integra l 
The matrix element is given by 


d 




y(x, - 


4it[(x - x) 
A ^ 

j 



x) 
2 + 



dx dy 


(34) 


where J, * 1,2, . . . ,M and j “ 1,2, . . . ,n 
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not Identical 


Coordinates of 1’, and P. 

i- 1 


If the coordinates of P,(x^,Vj^) and P.(x.,y.) are not identical, 

then equation (J4) involves a non-singular integral which can be evaluated 
at once. Substitution of the limits of integration yields 


yj+r. Xj+V 

/ / 

yj-qi Xj-6j 


y(Xj - x) 


[iXj, - x)2 + y‘ 


dx dy 


(B6) 


Hence 


where 


d|^j » -(1/16 ti) £n (N/D) 


y j + r .)2 + (x. - - (S.)2j _ q j2 + - 


(B7) 


D . [(yj - q.)- f ( 


N. * *j) 


X, - Xj - «^) 




+ r^ + (x^ - Xj, + 6, 


1 


are defined as shown in figure 1. 


The quantities 6 j , rj and qj 
Let 2hj be the height of A j , then 

2hj = rj + qj 

The quantities rj and qj are chosen as follows: 

r, . qj . hj 

if Pj is at the center of Aj ; or 

rj = 2bj and = 0 

if Pj is at the bottom edge of Aj; or 

rj “ 0 and qj = 2hj 
if Pj is at the top edge of Aj . 


(B8) 


(B9) 


(BIO) 


(BID 


Coordinates of P^ and Pj Identical 


The coordinates of P£(x£,yi^) and Pj(xj,yj) can be identical even 
if This occurs if P£ and Pj are both on the airfoil, ’lave the 

same x-coordinate , but are situated on opposite sides of the airfoil. 

If the coordinates of P£ and Pj are identical, equation (34) contains 
a single Integral. Following the method by Ogana (refs. 6 and 9), the 
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w 


Integral is evaluated after enclosing the singularity in a vanishing 
retangular cavity, (fig. 1); i.e.. 




y(Xj - x) 

- x)2 + 



dy 


(B12) 


Substitution of the lii.iits of integration and subsequent evaluation 
yields 



(B13) 


3. Matrix Associated with the Nonlinear Contribution to the Potential 


The matrix associated with the contribution of the nonlinear effects 
to the potential at a mesh point Pi(xj.yj^) is 




H - ^ 

4n - x^;2 + ( i; _ y^) J 


d(, dc 


(35) 


Diag o nal Elements 

If Aj contains then Pj = Pj and diagonal elements are 

obtained. '^The integral in equation (35) becomes singular at Pj^(x^,y£). 
Following the method by Ogana (refs. 6 and 9), integration is performed 
after enclosing the singularity in a vanishing rectangular cavity, Oj. 
(fig. 1). Equation (35) is consequently defined as follows: 



The appropriate integration steps yield 


bii - 0 


(B14) 


(B15) 


Off-diagonal Elements 

If Pj^ does not belong in Aj , off-diagonal elements are obtained 
by evaluating equation (35) written in the form 


17 


XI - f 


(IU6) 


Vi+r, Xj+6j 

-■ / / -r — T- 1 

/a .f. - x^)^ + (f. - yj) 

Vj~q^ Xj-6j L '• ^ J 


dr. 


Since* there are no sinRularit les , this Integral tan be evaluated 
at I’nce to Rive 

bij - -d/4n)[^Aj + Aj + A 3 + A^J, I i j 


wliere 


A, = 2 ^y.i - yi '■j> 


[ (Xj - + (Sj)‘ + (vj - Vj + i j) 

<:Xj - Xj - 5j)2 + (yj - Vj + rj) _ 


= 2 - Vi - qj) £n 


(Xj - xj - 6j)‘ + (yj - Xi - qj) 

(Xj - Xj + 5 j)- + (yj - Xj - qj)J 


(Xj - Xj + t'Sj) 


iro t; 


/•j - Yi + r,\ Aj - V, - 

anl : — T I ■ tanl 

V-j - -i + / V-'j - **1 


A4 = (Xj - Xj - iSj 


arc tan 


/Vj - Xj - qA /^ i - Vj + 

I 1 _ ;,i-c t..,nl 

\\1 - '‘I * ^:/ \ 1 ■ 


where rj and qj are chosen according to equations (K8) - (Bll) 


(B17) 
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Figure 2 
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AIRFOIL 


- Distribution of mesh points and rectangular elements in the 
upper half plane. 
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Figure 
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A.- Coefficient of pressure for a parabolic-arc airfoil of thickness 

ratio 0.06, at M “ 0.87. 
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